#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#ifdef __GNUC__
#include <iostream>
#include <fstream>
using std::ofstream;
using std::cerr;
#else
#include <iostream.h>
#include <fstream.h>
#endif
#include <string.h>
#if !defined(_Windows) && !defined(__linux__)
#include <graphics.h>
#include <conio.h>
#endif
#include "remez.h"

#ifndef __GNUC__
unsigned _stklen = 16384;
unsigned _heaplen = 0;
#endif

static char buffer[256];
static int len, bands, i;
static FilterTypes FilterType;
static SearchTypes SearchType;
static Filter f;
static FILE *input;
static ofstream output;
static long double a;
static Vector H;

#if !defined(_Windows) && !defined(__linux__)
static double MinY = -100, MaxY = 5;

static double WorldX1, WorldY1, WorldX2, WorldY2;
static int ScrX1, ScrY1, ScrX2, ScrY2;

void SetWorld(double x1, double y1, double x2, double y2)
{
  WorldX1 = x1;
  WorldY1 = y1;
  WorldX2 = x2;
  WorldY2 = y2;
}

void SetScreen(int x1, int y1, int x2, int y2)
{
  ScrX1 = x1;
  ScrY1 = y1;
  ScrX2 = x2;
  ScrY2 = y2;
}

int World2ScrX(double x)
{
  return ((x - WorldX1) * (ScrX2 - ScrX1) / (WorldX2 - WorldX1) + ScrX1);
}

int World2ScrY(double y)
{
  return ((y - WorldY1) * (ScrY1 - ScrY2) / (WorldY2 - WorldY1) + ScrY2);
}

double Scr2WorldX(int x)
{
  return ((x - ScrX1) * (WorldX2 - WorldX1) / (ScrX2 - ScrX1) + WorldX1);
}

double Scr2WorldY(int y)
{
  return ((y - ScrY1) * (WorldY1 - WorldY2) / (ScrY2 - ScrY1) + WorldY2);
}

static inline Float near CplxAbs(Float re, Float im)
{
  return (sqrt(re * re + im * im));
}

Float MFIR(const Vector &h, Float x)
{
  assert(x >= 0 && x <= M_PI);
  Float c2 = 2 * cos(x), x1, x2, x3;
  int k;
  for(k = Lo(h), x1 = x2 = 0; k <= Hi(h); k++)
  {
    x3 = h[k] - x1 + c2 * x2;
    x1 = x2;
    x2 = x3;
  }
  return (CplxAbs(x2 - 0.5 * c2 * x1, sin(x) * x1));
}

inline Float lnMFIR(const Vector &h, Float x)
{
  Float a = MFIR(h, x);
  return (Eq(a, 0) ? -1e10 : 20 * log10(a));
}

void ShowAmp(void)
{
  viewporttype vp;
  int driver = 0, mode = 0;
  const char *xaxis[] = { "0", "pi / 4", "pi / 2", "3 * pi / 4", "pi" };
  initgraph(&driver, &mode, "");
  if(graphresult() != grOk)
    return;
  getviewsettings(&vp);
  if(FilterType != OddDifferentiator)
  {
    vp.left += 40;
    vp.right -= 20;
    vp.top += 20;
    vp.bottom -= 30;
    SetScreen(vp.left, vp.top, vp.right, vp.bottom);
    SetWorld(0, MinY, M_PI, MaxY);
    setcolor(WHITE);
    settextstyle(DEFAULT_FONT, HORIZ_DIR, 1);
    settextjustify(CENTER_TEXT, TOP_TEXT);
    for(i = 0; i <= 4; i++)
    {
      moveto(World2ScrX(M_PI / 4 * i), World2ScrY(MinY));
      lineto(World2ScrX(M_PI / 4 * i), World2ScrY(MaxY));
      outtextxy(World2ScrX(M_PI / 4 * i), vp.bottom + 5, xaxis[i]);
    }
    settextjustify(RIGHT_TEXT, CENTER_TEXT);
    moveto(World2ScrX(0), World2ScrY(MinY));
    lineto(World2ScrX(M_PI), World2ScrY(MinY));
    sprintf(buffer, "%lg", (double)MinY);
    outtextxy(vp.left - 5, World2ScrY(MinY), buffer);
    for(i = floor(MinY / 10) + 1; i <= floor(MaxY / 10); i++)
    {
      moveto(World2ScrX(0), World2ScrY(i * 10));
      lineto(World2ScrX(M_PI), World2ScrY(i * 10));
      sprintf(buffer, "%i", i * 10);
      outtextxy(vp.left - 5, World2ScrY(i * 10), buffer);
    }
    moveto(World2ScrX(0), World2ScrY(MaxY));
    lineto(World2ScrX(M_PI), World2ScrY(MaxY));
    sprintf(buffer, "%lg", (double)MaxY);
    outtextxy(vp.left - 5, World2ScrY(MaxY), buffer);
    setviewport(vp.left, vp.top, ++vp.right, ++vp.bottom, 1);
    vp.right -= vp.left;
    vp.bottom -= vp.top;
    vp.top = vp.left = 0;
    SetScreen(vp.left, vp.top, vp.right, vp.bottom);
    setcolor(LIGHTRED);
    setlinestyle(SOLID_LINE, 0, THICK_WIDTH);
    moveto(World2ScrX(0), World2ScrY(lnMFIR(H, 0)));
    for(i = vp.left + 1; i <= vp.right; i++)
      lineto(i, World2ScrY(lnMFIR(H, Scr2WorldX(i))));
    setcolor(LIGHTGREEN);
    for(i = 1; i <= GetBandsNo(f); i++)
    {
      double y;
      y = Eq(f[i].GetTrans(), 0) ? MinY : 20 * log10(f[i].GetTrans());
      moveto(World2ScrX(f[i].GetLow()), World2ScrY(y));
      lineto(World2ScrX(f[i].GetHigh()), World2ScrY(y));
    }
  }
  else
  {
    vp.left += 90;
    vp.right -= 20;
    vp.top += 20;
    vp.bottom -= 30;
    SetScreen(vp.left, vp.top, vp.right, vp.bottom);
    SetWorld(0, 0, M_PI, M_PI);
    setcolor(WHITE);
    settextstyle(DEFAULT_FONT, HORIZ_DIR, 1);
    settextjustify(CENTER_TEXT, TOP_TEXT);
    for(i = 0; i <= 4; i++)
    {
      moveto(World2ScrX(M_PI / 4 * i), World2ScrY(0));
      lineto(World2ScrX(M_PI / 4 * i), World2ScrY(M_PI));
      outtextxy(World2ScrX(M_PI / 4 * i), vp.bottom + 5, xaxis[i]);
    }
    settextjustify(RIGHT_TEXT, CENTER_TEXT);
    for(i = 0; i <= 4; i++)
    {
      moveto(World2ScrX(0), World2ScrY(M_PI / 4 * i));
      lineto(World2ScrX(M_PI), World2ScrY(M_PI / 4 * i));
      outtextxy(vp.left - 5, World2ScrY(M_PI / 4 * i), xaxis[i]);
    }
    setviewport(vp.left, vp.top, ++vp.right, ++vp.bottom, 1);
    vp.right -= vp.left;
    vp.bottom -= vp.top;
    vp.top = vp.left = 0;
    SetScreen(vp.left, vp.top, vp.right, vp.bottom);
    setcolor(LIGHTRED);
    setlinestyle(SOLID_LINE, 0, THICK_WIDTH);
    moveto(World2ScrX(0), World2ScrY(MFIR(H, 0)));
    for(i = vp.left + 1; i <= vp.right; i++)
      lineto(i, World2ScrY(MFIR(H, Scr2WorldX(i))));
    setcolor(LIGHTGREEN);
    moveto(World2ScrX(0), World2ScrY(0));
    lineto(World2ScrX(f[1].GetHigh()), World2ScrY(f[1].GetHigh()));
  }
  getch();
  closegraph();
}
#endif //!defined(_Windows) && !defined(__linux__)



void InputError(void)
{
  cerr << "Error in input data\n";
  exit(3);
}

int main(int argc, char *argv[])
{
  if(argc >= 2)
    input = fopen(argv[1], "rt");
  else
    input = stdin;
  if(argc >= 3)
    output.open(argv[2]);
  else
#ifdef __linux__
    output.open("/dev/stdout");
#else
    output.open("CON");
#endif
  if(fscanf(input, "Filter length = %d ", &len) != 1)
    InputError();
  if(fscanf(input, "Filter type = %s ", buffer) != 1)
    InputError();
  if(strcmp(buffer, "OddSymetrical") == 0)
    FilterType = OddSymetrical;
  else
    if(strcmp(buffer, "OddAntysymetrical") == 0)
      FilterType = OddAntysymetrical;
    else
      if(strcmp(buffer, "EvenSymetrical") == 0)
        FilterType = EvenSymetrical;
      else
        if(strcmp(buffer, "EvenAntysymetrical") == 0)
          FilterType = EvenAntysymetrical;
        else
          if(strcmp(buffer, "OddDifferentiator") == 0)
            FilterType = OddDifferentiator;
          else
            InputError();
  fscanf(input, "Search type = %s ", buffer);
  if(strcmp(buffer, "ExhaustiveOnly") == 0)
    SearchType = ExhaustiveOnly;
  else
    if(strcmp(buffer, "UseSelective") == 0)
      SearchType = UseSelective;
    else
      if(strcmp(buffer, "UseSelectCubic") == 0)
        SearchType = UseSelectCubic;
      else
        InputError();
  if(fscanf(input, "Bands = %d ", &bands) != 1)
    InputError();
  f.SetNoBands(bands);
  for(i = 1; i <= bands; i++)
  {
    sprintf(buffer, "Lower edge %d = %%Lg ", i);
    if(fscanf(input, buffer, &a) != 1)
      InputError();
    f[i].SetLow(a);
    sprintf(buffer, "Upper edge %d = %%Lg ", i);
    if(fscanf(input, buffer, &a) != 1)
      InputError();
    f[i].SetHigh(a);
    sprintf(buffer, "Transmitance %d = %%Lg ", i);
    if(fscanf(input, buffer, &a) != 1)
      InputError();
    f[i].SetTrans(a);
    sprintf(buffer, "Error %d = %%Lg ", i);
    if(fscanf(input, buffer, &a) != 1)
      InputError();
    f[i].SetErrorRatio(a);
  }
  if(fscanf(input, "Grid = %u ", &S) != 1)
    InputError();
  if(fscanf(input, "Required Q = %Lg ", &a) != 1)
    InputError();
  RequiredQ = a;
  if(argc >= 2)
    fclose(input);
  InitRemez(f, len, FilterType);
  if(Remez(SearchType, 20) != 0)
    cerr << "No convergence\n";
  else
  {
    H = CalcH();
    for(i = Lo(H); i <= Hi(H); i++)
      output << "h[" << i << "]\t= " << H[i] << "\n";
#if !defined(_Windows) && !defined(__linux__)
    ShowAmp();
#endif
  }
  return (0);
}
